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We use a path integral approach to calculate the superfluid density of a Bose lattice gas in the 
limit where the number of atoms per site is large. Our analytical expressions agree with numerical 
results on small systems for low temperatures and relatively weak interactions. We also calculate 
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develop tools for calculating discrete time path integrals. These tools should be broadly applicable 
to a range of systems which are naturally described within an overcomplete basis. 
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I. INTRODUCTION 

Superfluidity is one of the most profound collective 
manifestations of quantum mechanics [TJ [5] . It is char- 
acterized by dissipation-less flow and is analogous to the 
vanishing resistivity seen in superconductors. The phe- 
nomenology of superfluidity is largely contained in Lan- 
dau's two fluid model: one component, the normal fluid, 
responds to the motion of the container walls, while the 
other component, the superfluid, does not. The total 
density p = p n + p s is the sum of the density of each 
component. Leggett showed that at zero temperature, in 
a translationally invariant system, either p s = or p n = 
[3]. In a lattice, however, even at T = 0, p s /p n can be 
finite. Here we calculate the superfluid fraction for an 
interacting Bose lattice gas in the large filling limit. Our 
study complements continuum calculations of superfluid 
densities [IHZ]. 

We are largely motivated by experiments of cold 
bosonic atoms in optical lattices j8]. These systems are 
well described by the Bose-Hubbard model jS] , which can 
be studied using mean field theories [10] and Quantum 
Monte Carlo methods [TTJ U2] . Further motivated by ex- 
periments where two bosonic species are trapped on a 
lattice [T3"HTrj] , we also calculate the superfluid density of 
a two-component system. Such mixtures have rich be- 
havior, including exotic phases such as paired superflow 
and counter-superflow |16[ I17j . 

To calculate the superfluid fraction we use a func- 
tional integral approach where we include quadratic fluc- 
tuations about a coherent state which makes the action 
stationary. This method becomes exact in the weakly- 
interacting, low-temperature, high-density limit. We give 
finite temperature results and compare with exact nu- 
merical diagonalization on small systems. 

Our calculation involves coherent state path integrals. 
As was previously established [T5] there are difficulties 
with the continuous time limit of these objects. We show 
explicitly how to calculate the discrete time path inte- 
grals. The resulting formalism contains extra terms not 
seen in the standard approach. 

In Section [TT] we introduce the physical meaning and 
thermodynamic definition of the superfluid density. In 
Section |III| we present the results of the calculation in 



the case of a single species of bosons on a lattice, and in 
Section [TV we explore the superfluid properties of two- 



component bosons. Appendix [A] highlights the necessity 
of the discrete-time formalism in the use of coherent state 
path integrals, and Appendix [B] demonstrates the tech- 
nical details of calculating thermodynamic quantities in 
this formalism. 



II. SUPERFLUID DENSITY 

To define the superfluid density p s we follow |2 and in- 
troduce a new thermodynamic variable v s via a thought 
experiment. We imagine a fluid at rest within an in- 
finitely long cylinder that is itself at rest. This defines 
the lab frame. We now give the cylinder an infinitesi- 
mal velocity — v s along its axis. After we have allowed 
the container and fluid to reach equilibrium, the mass 
current as observed in the cylinder frame of reference is 



3 = PsV s 



(1) 



which defines p s , the superfluid density. A normal fluid 
will move as a rigid body with the container and so have 
p s = 0; an entirely superfluid liquid will feel no drag and 
remain at rest in the lab frame, yielding p s = p. It is also 
convenient to define the normal density, 



P = Ps + Pn- 



(2) 



Formally, we may calculate the superfluid density as 
the second derivative of the free energy density T with 
respect to v s , 



d 2 T 



dvl 



(3) 



v s =0 



In a more technical language, this indicates that the su- 
perfluid density is the low- frequency, long wavelength 
limit of a transverse current-current correlation function 

In a translationally invariant system, for a fluid with 
well-defined quasiparticles, one can express Eq. (J3| as a 



2 



sum over the excitation spectrum, 



Hamiltonian, 



Pn = 



d 3 p 

(2ttK) 2 



P ■ v 6 



dn b 



(4) 



v.=0 



where rib — \eP £k — l] is the Bose-Einstein distribu- 
tion function and e p is the energy of an excitation of 
momentum p. 

In three dimensions, the microscopic understand- 
ing of superfluidity involves condensation into a single 
macroscopically-occupied quantum state. If the wave- 
function of that condensed state is given by ip(r,t) — 
\J ' p c (r,t)e tx ( r,t ' , then the superfluid velocity v s is di- 
rectly related to the phase x, 



-Vx(r,t). 



(5) 



The variable p c defines the condensate fraction, p c j p, the 
portion of the system that is condensed into the ground 
state. This fraction is not, in general, equal to the super- 
fluid fraction p s /p. 



Experimental probes of p 3 

To measure p s in a gas of cold atoms we propose the 
following experiment. One begins with an equilibrated 
Bose gas in an optical lattice, confined by an additional 
harmonic trap. The dimensionality can be controlled by 
adjusting the intensity of the lattice beams in the rele- 
vant directions. The harmonic trap is then turned off, 
and the lattice accelerated to velocity v s by chirping the 
frequency of one of the lattice beams. One then turns off 
the lattice and uses time-of-flight expansion to measure 
the momentum p of the cloud. In the limit that all steps 
are adiabatic, the mass contained in the normal compo- 
nent is p/v s . Converting this to a density or a superfluid 
fraction is trivial. 

Gadway et al [3T] have implemented a related protocol, 
but did not emphasize the fact that they were measuring 
the superfluid density. Alternate theoretical proposals in- 
volve rotation or artificial gauge fields. Ho and Zhou [22] 
showed that the superfluid density can be extracted from 
images of rotating clouds. John, Hadzibabic and Cooper 
[2"5] identified a global spectroscopic measure of super- 
fluidity, while Carusotto and Castin [24] investigated a 
local probe. 



III. SINGLE SPECIES 
A. Model 

We begin by analyzing the case of a single species of 
weakly-interacting bosons on an optical lattice. Such as 
system can be modeled by the single-band Bose- Hubbard 



H = - 



(id) 



E 



u 



(hi - 1) - phi 



(6) 



where summations are over the sites i and over the pairs 
of nearest neighbors (i, j) . Here a, (aj) is the annihilation 
(creation) operator for a boson on site i and hi = djdj is 
the number operator for the site. In this paper we focus 
on the case of a cubic D-dimensional lattice, taking the 
lattice spacing to be ao and the volume of the system to 
be V. 

The Bose Hubbard model is a good description of 
the atomic system as long as the band spacing Eb is 
greater than all relevant energy scales in the system, 
Eb 3> J, U, T. Under these conditions, excitation into 
higher bands can be neglected. In cold atom exper- 
iments this spacing scales as Eb ~ ^/£VqEr where 

h 2 h 2 

Er — is the recoil energy for particles of mass 

m trapped by lasers of wavenumber k — 2n/\, and Vq 
is the optical lattice depth, which is typically of order 
V ~ 10— 100 x Er. For near-optical lasers and particles 
lighter than m < 100 amu the single band approximation 
works up to T < 10~ 6 if [25] . 

We introduce the velocity v s into our model by apply- 
ing a phase twist AO to the hopping term, 



a,- o,- — > e 



-iAO.(r 1 -r 1 )/aa a J 6 



(7) 



or equivalently, <L- — > e lAL@ ' rj ^ a °a J where r, is the posi- 
tion of lattice site i. This phase is related to the lattice 
velocity by 



-A© 



mao 

and so wc obtain the relation 



dd' 



n 2„2 



n 2 



d 2 T 



dAQ d dAe d , 



(8) 



(9) 



A©=0 



where d, dl = 1, . . . , D = x, y, z are the lattice directions. 
In principle, the superfluid density on a lattice may be 
a symmetric rank 2 tensor, but for the cubic lattice one 
has p d s d ' = Sdd'Ps- 

Like all thermodynamic quantities, the free energy 
density can be derived from the partition function, 



given by 



Z = Tre-^=5>|, 

I*) 



-0H 



(10) 



(11) 



where (3 = 1/T is the inverse temperature and the 
sum is over a complete set of states Introduc- 
ing the overcomplete coherent state basis, a,i\pi,ipi) = 



3 



y/fhe lipi \pi,(fi), we break up the operator e~^ H into N t 
slices and express the partition function as a path integral 
of the Euclidean action over the classical fields l26l. 



Z = <j> DpDip exp [— Se] ■ 



(12) 



As discussed in Appendix [A] one must use the discrete 
time formulation of the action, 



JV t -l 

t=o 



(13) 



with 



l e = ^2 - lo S [(Pi,t, <Pi,t I Pi,t+l,fi,t+l)} 
i 

N t (pi,t,fi,t | Pi,t+i,fi,t+i) 

= £ PM \ Pt ' t+1 - Vp^w^ t+1 -^ t] 



JAt ^2VPhtP],t+ie 



i(tpj.t+i-<Pi,t-A6ji) 



(14) 



+ \/PjJPi~t+ie 



i(<Pi,t+l — <Pj,t— A0y) 



UAt 



y^,Pi,rPi,T+ie 



-<Pi,r) 



pAt V/VPi,r+ie i( ^ T+1 - Vi > T) 



where A0y = A© • (r, — r\,) /ao and At = f3/N t is the 
discrete time step. We take the number of time steps Nt 
to be large. 



B. Saddle-point Approximation 

We expand the fields pi, (pi around the mean density 
p and mean phase twist A<I> = ^A^f,^, with the 
unit vector in direction d. For any site i and its nearest 
neighbors along d, i+d and i-d, we have 



Pi,t — P + Spi,t 

'Pl.t = —ri ■ A# + <f> i t 
a 

<Pi +d ,t - <Pi,t = A$ rf + 4>i +d , t - <t>i,t 

<p%,t - <Pi- d ,t = A<i>d + 4>i, t - i>i_ t ,f 



(15) 



We take these perturbations to be small, Spij <C p, 
<t>i± d ,t - 4>i,t < 1) 4>i,t+x - (t>i,t < 1. The v alidity of 
these assumptions is examined below, in Sec. |IIIE| In 
particular, when T,U < pJ one finds (5p? t \ ~ p and 

(^(4>x+i — </>:r) 2 ^ < 1//0- Thus if p <C 1 this expansion is 
well behaved. 

Although we assume ((pi ±d ,t ~ <t>i,t) and ((f> it t+i - <f>i,t) 
are small, we make no assumption that (pa itself is small. 



Consequently our calculation is valid even in low dimen- 
sions, where the condensate fraction vanishes and there 
is no long range order. 



Eq. ( 14 1 expanded around the mean values reads 



^E 



(16) 



where each subsequent term involves higher powers of the 
fluctuations. 

The first term is a constant, 



Y,d 2 P Jcos (A$ d - A9 d ) 



+f f - PP 



At. 



(17) 



Keeping only this term gives the mean-field Gross- 
Pitaevskii approximation where p s = p = p. 
The second term, linear in the perturbation, is 



-2J£ d cos(A$ d - AO d ) 
+Up- p 



Atdp l 



(18) 



The saddle-point mean values minimizing Cq are 



A$ = A6 



JJ ^ + 2j£cos(A$ d )^J 



(19) 



Setting p to this value makes C\ vanish. Such a structure 
is generic, as minimizing the zeroth-order action causes 
the first order action to vanish. To calculate the super- 
fluid density, we take A© = but keep A3? finite, giving 
the bosons velocity — — A<I> relative to the lattice. The 

superfluid density becomes p s = 

n |_ozi<P d j ^,£_o 

The "interaction" term, which we neglect in our calcu- 
lations, consists of terms of third order or higher in the 
perturbation fields, 



£\nt/P = 0(8p/p,<t>j - (pi) 



3 = 0(1/ VP) 3 - 



(20) 



Our non-trivial results come from the the quadratic 
term, which is best expressed in momentum space, 



a?d D k 



(27T) 1 



V 



«o 



D 



k,uj n 



(21) 



_ N t -1 
2 



with fre- 



where summation is over n 
quencies given by ui n = n^n, and the integration is over 
the first Brillouin zone \kd\ < Tt/a . 

Appendix [B] provides details on the explicit form of 
C 2 ,u and the calculation of propagators. However, all 
significant physical results rely only on the behavior of 
the propagators and action at two regimes: ujAt -C 1 
(superscript p for pole behavior) and uAt = ire lx (su- 
perscript o for contour behavior). These are given, at 
A$ = 0, by 



4 



V 



( 5 P 5 P)i,u> = ~dP 

Lip. 



1 2E\k 



(<W> p fc ,. = 



v 



At 


io 2 + £ k 2 


1 




At 


io 2 + £ k 2 


1 


2^2fc 



V 1 

afiip [Atu 2 +e k 



+ o {Aty 

+ O {At)° 
+ O {At)° 



V 



( 5 P s P)l, u = ~dP 



<<W>)fe, w = n 

LI. 



V 



o 

V 1 

a? 4p 



1 
1 

21 
1 + 



£ lk At 



1 — cos (7re lx ) 

sin (-7re* x ) 
— cos (7re J x) 

1 — cos (7re lx ) 



0(At) J 



o {Aty 



o (Aty 



(22) 



where we use the notation (XY) k u = (Xk, u Y-k,-u) an d 



on the right it is understood u) = 
appearing in these expressions are 



£ik =4J^sin 2 (k d a /2) 

d 



The energies 



-2k 



2pU + ( 4J^sin 2 (k d a /2) 



P 2 — 



4J^sin 2 {k d a /2) 

d 

\ " 



(23) 



2pU+ 4J> sin^ (k d a a /2) 



d 



In the continuum limit, one has £] k — > ^V 5 - and £ 



-2k: 



i 2gp, where g — Ua^ /m, p = m(n) faff. The 
excitation spectrum £ k then corresponds to the familiar 
Bogoliubov result. 



C. Superfluid Density 

The superfluid density is given by 



n 2 





~d 2 lnZ~ 




(-w) 


[dA<t> 2 _ 


A#=0 



(24) 



Some insight may be gained by inserting the path integral 
expressions for the free energy density and the partition 
function into this equation. We find that the superfluid 
density, to order O (1/p) , can be decomposed into three 
terms, 



Ps 



2maQj m 



r, 2 



[no 



(25) 



We identify the first term as the total density, from which 
two normal-density terms are subtracted. 
These are, respectively, 



2J 



1 N t V 



d 2 £ 



dA<5> 2 



-T 



a?d D k 



d 2 p dC k 2 



(27T) 



dA<P 2 



dp 



(26) 



2J v 

„D U n 



1 

W 



E 





1 d 2 c k 2 u - 


/ (2,f < 


\ dA^> 2 



(27) 





m,n 




dCf Wm 


\ SA* d 


_ dA4> d 
p 



a"d D k a"d D q 
(2nf (2tt) D 



0A$ rf 



(28) 



where (X) = ~ § "DpVipX exp[iS#]; on the right-hand 
side of the equations, derivatives in A§ d and p are to 
be taken at a constant p and A<fr, respectively, and then 
evaluated at A<& = 0; and we have omitted multiple 
vanishing terms. This is similar to the calculation shown 
explicitly in Appendix [B| 

The term tiq = (n) = — |f is the total average occu- 
pation number. It is given by 



n = p + 



1 



a D d D k 
(2tt) D 



£lk 

~£~k~ 



coth(/3£ fc /2) . (29) 



At T = 0, one finds n -> (/i + 2D J) /U as J7/J 0, 
and n — > /i/f/ + ^ as C//J — > oo. These correspond 
to the correct occupation numbers in the non-interacting 
and the no-hopping regimes. The explicit calculation of 
this term is given in Appendix [B] 
The term nf^ is given by 



n 



Sp<f> _ 



a u d u k 
(2nf 



2Jsin^ (k d a ) - 



dn b \ 

d£ k J A<s>=0 



(30) 



This expression is reminiscent of the form of the normal 
density in the continuum case, given in Eq. Q, and it 
likewise vanishes at T = 0. 



The additional term, 



can be understood to come 



from density-density and phase-phase correlations cre- 
ated by the interaction term in the hamiltonian. 



a D d D k 

=- (1 - cos (k d a )) x 

(2tt) D 



{£lk + £lk) 

2£ k 



coth (/3E k /2) - 1 



(31) 



A*/ 



At T = 0, U = 0, this term vanishes and the superfluid 
fraction becomes one; at non-zero values of U this term 
is finite even at T = 0. 



5 



The resulting superfluid fraction, p s /p, is plotted in 
Fig. JT] at zero temperature as a function of U /J and in 
Fig. pffor set values of U / J as a function of temperature. 
In Fig. [3] we show the curve U (T) where p s vanishes, 
suggesting a phase transition. The limits of validity of 
these results will be discussed in Sec. IIIIEI 




u/j 



FIG. 1: (Color online) The superfluid fraction p a /p as func- 
tion of U/J in an infinite 3D cubic lattice, for (n) = 10 (solid 
red line) and (n) = 1 (dashed blue line), calculated to lead- 
ing order in a 1/ (n) expansion. As discussed in the text, the 
results are not expected to be quantitatively accurate above 
U/J>{n). 



pjp 

1.0 c 



FIG. 3: The values of U/J,T/J at the intercept p a = 0, 
suggesting a superfluid-Mott insulator transition. The cal- 
culation is performed for an infinite 3D cubic lattice, with 
(n) = 10. The inset shows the form of the curve at small 
T/J. As discussed in the text, the results are not expected 
to be quantitatively accurate at values of U/J > (n) or 
T/y/ J (J + pU) > (n), but the form is qualitatively similar 
to curves generated by quantum Monte Carlo methods [27] . 



First we compare our result to the continuum limit 
by taking ao — > 0, J — > 00 so that Joq is constant. In 



this case, the second term Eq. (31) vanishes. This can 



be seen by separately considering the contributions from 
k <~*j 1/ao and k <C 1/ao; in both cases the integrand 
vanishes. The first part of the normal density, Eq. (30), 



has contribution only from finite momenta, and becomes 



n£'* -> 2m Ja, 



n '"n 



d D k 2 f_dUb 



(27T) 1 



(32) 



A$=0 



T/J 



FIG. 2: (Color online) The superfluid fraction p a /p as func- 
tion of T/J in an infinite 3D cubic lattice, for (n) = 10, at 
U/J = 0.01 (dotted blue line), U/J = 1 (dashed red line) 
and U/J = 100 (solid yellow line). At U = 0, p s vanishes 
at T/J = 41.5, the ideal gas transition temperature. As dis- 
cussed in the text, the re sults are no t expected to be quanti- 
tatively accurate at T / y/ J (J + pU) > (n). 



D. Analytical Limits 



Here we examine the behavior of Eq. (25) in several 
limiting cases. 



with the continuum spectrum. Identifying 

2mJa^/h 2 — > 1, this is precisely the known continuum 
result seen in Eq. Q. 

Another important limit is zero temperature and large 
p. If U ~ J, then pU ^> J and to first order w 

yj 8U Jp (J^d sin2 (kd.ao/2)) , £ik + £2k ~ 2pU and the 
normal density becomes 



DjDi 




where ^= 



= / 



cos (k d a )) 
1 

~ 2 

(l-COSfgg)) 



1 



1 



(33) 



1 

/20' ' 



1 

/35 : 



1 

'51 



( 2 ^) D ^32(^,3^(6^2)) 

in one, two and three dimensions respectively. This ex- 
pression is suggestive of a phase transition from super- 
fluid to Mott insulator at U — U c ~ fnpJ ■ In two and 
three dimensions, the values for fjj are about double the 
mean- field result of U c - 2D x 4 (n + ±) [28j- More 
comparisons along these lines are made in Sec. |III F| 
As discussed in Section III E these estimates are beyond 



6 



the range of U/J where our approximations are quanti- 
tatively valid. It is nonetheless appealing to see the Mott 
transition appearing within this formalism. 

Finally we consider the free particle case U = 0. There 
we have as a function of T 



where the first constraint is universally required, the sec- 
ond stems from the density fluctuations in Eq. (381 and 



the last two from the phase fluctuations in Eq. ( 37 1 



F. Gutzwiller Ansatz 



+ 



a D d D k 1 



(2tt) D 2 

D d D k 1 
2 



ao 



(27T) 



D 



(coth(/3£ fe /2)-l) 
cos (k d a ) (1 - coth (/3£ k /2)) 

. fij sin 2 (fc d q ) 
S inh 2 (/3£ fc /2) 



(34) 



The integrand in the first line \ (coth (/3£fc/2) — 1) = 
rib (£k) and the one in the second is a total derivative 
that vanishes at fc^ao = 0, 2n, and so n n = n ex , the total 
occupation of excited states. p s vanishes at (n) — n ex , 
corresponding to the ideal gas transition temperature. 



At zero temperature, an alternative approach to cal- 
culating the superffuid density is to use the Gutzwiller 
ansatz [291 . This is an uncontrolled variational method 
which reproduces the Bogoliubov results at weak cou- 
pling |10j . and gives us a point of comparison for our 
results. 

The Gutzwiller approach assumes that the ground 
state of the lattice system may be decomposed into a 
product of single-site states, 



(40) 



E. Realm of Validity 



where 



Though the formulation of the action in Eq. (14), (21) 



is exact, our calculations are performed by neglecting the 
infinite series of terms in C^. We can place bounds on 
the realms of validity of this approximation by requiring 
that the perturbations from the mean values p, A$ be 
small, 



{5pi,t$Pi,t) ^ P 2 
k +d ,t - <kt) 2 



< 1 



(35) 



i,t+l 



We do not require the phases themselves to be small, only 
the deviation from one site to another and from one time 
step to another. 

These fluctuations can be calculated by the use of the 
propagators in Eq. (22), 



k,f) = ~ 



l 



^^coth(/?4/2) 

(27T) J 



(SpiSpi) 



1 



,DjD 



d D k £u 

£k 



(2n) 1 



cabh(J3£ k /2) 



(36) 



(37) 



(38) 



Examination of the integrals in the latter two inequalities 
implies that to to keep these parameters small we must 
have 



P>1 
U/J < p 
T/J<P 



(39) 



T/^/j(J + pU)<p, 



\g) t 



oc 

£< 

m=0 



(41) 



and |0)j denotes the zero-boson state of site i. One finds 
the coefficients a™,^ m by minimizing the expectation 
value (i(jg\ H \ipo) ■ We expect a™ = a m to be equal 
on all sites and 9™ = (A3? ■ ri/ao)m where A$ is the 
phase twist, as before. 
One then minimizes 

(i/>g\B\1>g) = 

^2 - 2J^]cos (A$ d ) x 



J2 Vr^+la m+1 a m 



(42) 



+ 



' u < 

— m (m 
. 2 ^ 



1) 



pm) |a m | 5 



and finds that the superfluid density is 
d 2 



Ps 



m 2 a 2 1 
K 2 V 



d^ 2 



u \ m 



2ma 2 J m 



(43) 



la m+1 a m 



We calculated the parameters a m numerically by cut- 
ting off the sum at rn = 20. We compare the results with 
those of Eqs. ((25j)-((3T]> in Fig. [2} 



G. Numerical Comparison 



We also compared the results of Eqs. (|25[)- ( 31 1 to an 



exact numerical calculation of the superfluid density for 
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FIG. 4: (Color online) The superfluid fraction p a j ' p for an in- 
finite three-dimensional cubic lattice with (n) — 10, at T — 0. 
The dashed blue line shows the result of the Gutziller-ansatz 
calculation and the solid red line shows result as calculated 



using Eqs. (25 h(31 1 



a variety of small lattices in one and two dimensions. 
For a finite lattice and fixed number of particles, we can 
represent the Hamiltonian in Eq. ^ as a finite matrix. 
We diagonalized this matrix, finding all eigenstates and 
eigenvalues. We calculated the superfluid density by per- 
forming the full weighted trace over all eigenstates. 

We find that at zero temperature the approximate an- 
alytic expressions for the superfluid density match the 
numerical result well even at a relatively small number of 
particles per site, (n) = 4. Moreover the agreement per- 
sists to relatively large U. One such example is shown in 
Fig. [5] The finite temperature values do not agree as well 
with the numerical result, except for very large values of 
(n), but they follow the same trend as the numerically 
calculated results (see Fig. [6|. Overall the numerical re- 
sults confirm the limits of validity in Eq. ( 39 ) . 



For these comparisones we replaced the integrals in 
Eq. ( 29 )-( 31 ) with sums, corresponding to the finite size 



system. 



FIG. 5: (Color online) The superfluid fraction p 3 /p for a two- 
dimensional two-by-two lattice with 16 particles, at T = 0. 
The dashed blue line shows the numerically exact result and 
the solid red line shows the analytic approximation. 
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IV. TWO-COMPONENT SYSTEMS 



FIG. 6: (Color online) The superfluid fraction p s /p as a func- 
tion of temperature, for a two-dimensional two-by-two lattice 
with 16 particles, at two values of U/J. The dashed blue 
line shows the numerically exact result and the solid red line 
shows the analytic approximation. 



A. Model 

We apply the same path integral method to a system 
of two species of bosons on a lattice. The Hamiltonian 
for this system is given by 



H — -\- -\- H 



n 



(44) 



where H a for a =t, I are the single-particle Hamilto- 
nians for particle species f, i respectively, identical to 
Eq. ([6]) except with constants J a ,U a ,p (7 and operators 
a?, (af)^ as appropriate. The final term 



(45) 



is the inter-species interaction Hamiltonian. 
The action for this Hamiltonian given by 



Se — 



E 



E 



(£) J 

V 



a 

, 2 



D 



d°k 



r£-0 + cr-^V 



(46) 



where aCp^.int are again identical to those defined in 



Eq. ( 17 ) , ( 20 ) , ( Bl I with mean densities and phase twists 



p a , substituted for the one-component equivalents 



8 



as appropriate, and the saddlepoint relation 
p a = U a p a + U n p s - 2J a cos (A#5) 



used, with a indicating the non-cr species, so that j" =\ 
, \ =\. The additional terms are 



(47) 



n C k '" = 2 x -U n At 



l±cos^At)\ X J c 4- 



8pL,8pt h ,- u - Wl (1 - cos (wAt)) 4 iU 0t 



+ sin (wAt) PiSpl^t^ + sin (wAt) P t H>-fe, 



.4- j.t 



(48) 



and the higher order terms scale as ^^C^ k = p^p^U^i x O (l/^p) 3 . 
The in-species propagators are now, at wAt <C 1, 

V 



2 , c . - ! 

At ( w a + f +fc 3 Hw 2 + f_ fc 



1 



O (At) L 



o 

V 1 



At (oj 2 +£ +k 2 ) (oj 2 +£_ k 2 ) 

" 1 2£ 2fc (a; 2 + £ gfc 2 ) - 8p T ^£/ n 2 g g i fc 
At ( W 2 +£ +fc 2 ) (w 2 +£_ fc 2 ) 



(49) 



O (At) c 



while the contour pieces are identical to the single-component case. £ a k-,£a\k-,£a\k are the single-particle dispersion 
relations given in Eq. (23), and the new dispersion relations are given by 



£±k 2 — 



± 



c 2 a 2x 2 



4ptPiU^i 2 £^ik£iik- 



(50) 



The interspecies propagators are given by 



V 



V 



"o 



1 




^PtPi£n £ n 




At 


(uj 2 


+ £+ 2 ) (u 2 + 


£J) 


1 




2p a £\ a u) 




At 


(iO 2 


+ £+ 2 ) (c 2 + 


£J) 


1 




UJ 2 




At 


(oj 2 


f £ + 2 ) fw 2 + 


£- 2 ) 



O(At) 



0(At)° 



<3(At) f 



V 



k,uj 



o(Aty 



v 



<W> fc =-^n o(At) 



(51) 



V 



1 sin 2 (fe 2 *) At 

2 (1 - cos (Tre^)) 2 



O(At)' 



B. Superfluid Density 

In the presence of two species there are now three su- 
perfluid densities, 



p7 = 



n 2 



d 2 T 
dA^ a d dA^ a 



(52) 



where p° T is the superfluid response of species a to the 
twisting of the phase of species r. The diagonal terms 
p" s a are the superfluid densities of species cr, while the 



off-diagonal term p 



p. 



is the cross-stiffness. 



The full expressions for all three terms may be calcu- 
lated in a similar manner to the single-species case, as 
described in Appendix [EJ At zero temperature, the su- 
perfluid densities are given by 

^ 2m a a\j a m a 



r, 2 



,aU1 



(53) 



where the number of normal atoms per site is given by 

D (1 — COS (fcjOo)) x 



1 f a£d D k 



(2nY 

{£<jik+£ e j2k)(£+k£-k-\-£<Tk 2 )-^p\PiU^\ 2 £ a 
2£+k£-k(£+k-\-£-k) 



(54) 



1 



9 



The cross-stiffness at T = is 



n n - 



agd 3 k 
(2tt) 3 



ft 2 



sin 2 (& d ao) 



£+£_(£++£-) 3 



(55) 



While the cross-stiffness is expected to be substan- 
tial in hard-core bosons [T7] it is negligible in the weak- 
interaction case, as illustrated in Fig. [7] 

A more dramatic effect can be seen in the superffuid 
densities. At zero temperature, a strong coupling to a 
second species of particles can replenish the superffuid 
fraction, as long as the superffuid has an equal or larger 
hopping parameter. This is seen in Fig. [8] 




FIG. 7: The superffuid cross-stiffness pIV p^ as function of the 
interspecies interaction Ufi/U-f for a two-component Bose gas 
on an infinite 3D cubic lattice. Here (n ) = (n ) = 10, with 
Ui — Uf = 10 Jj, = 10 Jf, calculated to leading order in a 
1/ (n) expansion. 



V. OUTLOOK 

The T — normal density, p„, for lattice bosons is 
generally non-zero. As discussed in Sec. |Hf[ this prop- 
erty, and the temperature dependence of the superffuid 
density, can be experimentally studied using cold atoms. 
Here we calculated p n , and proposed comparing our re- 
sults with experiment. 

For our calculation we extended the standard saddle- 
point functional integral approach. When using a coher- 
ent state basis, the discrete time path integral contains 
extra terms over the continuous time limit version. We 
explicitly derived those corrections for the Bose Hubbard 
model. Similar issues appear in spin models, and our 
techniques could be applied there. 

Our results are applicable at high density, low temper- 
ature, and weak interaction. One could envision extend- 
ing them to strong interaction by using a different set of 
coherent states. For example, in the hard core limit it 
would be natural to use \9,<p)j = cos 8 \0) i + e lip sin 




i — , — , — , — , — i — , — , — , — , — l — , — , — , — , — i — , — , — , — , — i Ut/Uu 
-1.0 -0.5 0.0 0.5 1.0 

FIG. 8: (Color online) The superfluid fraction of the one 
component pt/p as function of the interspecies interaction 
U-fi/Uf for a two-component Bose gas on an infinite 3D cubic 
lattice. Here (n r ) = (n l ) = 10, Ui = U t = 10J t , and the 
second component has hopping parameter Ji/Jf = 0.1 (dot- 
ted blue line), J^/J-f = 1 (dashed red line) and J±/Jf = 10 
(solid yellow line). 



where (0)^,1 1), are the states with no particles or one 
particle on site i respectively. The other approach to ex- 
tending the validity of our results would be to include 
perturbative corrections. In particular, one might envi- 
sion summing an infinite set of these corrections using 
Feynman diagram techniques. 

We also present results for the superfluid properties of 
two-component lattice bosons. These are an active area 
of research, and there are rich possibilities for exploring 
our formalism in those systems. One experiment [14] has 
seen hints of the impact of one bosonic species on the 
superffuid properties of another. Those results appear to 
be in the opposite direction from our predictions - how- 
ever they are in a stronger interacting regime, near the 
superfiuid-to-Mott insulator transition and the quanti- 
tative applicability of our results to their experiment is 
questionable. We also neglect any processes which in- 
volve higher bands. 
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Appendix A: Discrete Time Path Integrals 

The traditional path integral formulation of quantum 
mechanics involves the transformation of the partition 
function 



Z = Tre-^ = J2' 

\4>) 



-m 



(Al) 
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into a path integral. Here {\ip)} is any complete basis of 
the states. In our case we will use the overcomplete basis 
of coherent states. 

To make the transformation, we break up the operator 



-PH/Nt 



into Nt time steps. We then insert 



an identity operator 1 = X) 



Z= £ nWe-^hfefi) 
{|V>*» * 



between each step, 
(A2) 



where the summation is over N t copies of the the basis 
\t/j), the product is over t — O..N t — 1 and we define 
\lpN t ) = |V>o)- 

Taking the number of time steps Nt <C 1 to be very 
large, we then expand 



(ipt I t/Jt+i 



1 iV't+i) 
Nt 



^ t \H\^ + i)+0{P/Ntf 



(A3) 



exp 



log [{ipt | ipt+i)} 



P (i/j t \H\Tp t +i) 
N {ipt I it+i) 



The integration over all Nt values of \ipt) is a path inte- 
gral, and one finds 



(A4) 



Z = J ViP e 
where Se — X)t an d 
Lt = - log [{ipt | Vt+i 



f3 {ip t \H\ipt+i) 



N {tpt | ip-. 



(A5) 



t+i) 



It is here that we diverge from the tradition continuous 
formulation of the path integral. One typically assumes 



l + (P/Nt)dt + 0(P/N t ) 2 )\iPt), (A6) 



and by taking Nt — ► oo the sum can then be converted 
into an integral, S™ nt — J dtL t , where 



Lt = -{iP {t)\d t \iP{t)) + {iP {t)\H\iP{t)) 



(A7) 



The expression in Eq. ( A6 ) is not always valid. In par- 
ticular, for an overcomplete basis the overlap {ip t \ ipt+i) 
remains finite for states that differ to a non-infinitesimal 
degree, and the difference between \ipt+i) and \ipt) need 
not go to zero as Nt — > oo. This leads to a breakdown of 
the traditional continuous path integral, as shown in |18j . 



However, the discrete time formulation Eq. (A5| remains 
valid. 

One example is the single-site Bose-Hubbard model, 



H ss = (ft - 1) - //ft 
where n is the number operator. 



(A8) 



The partition function for this Hamiltonian can be cal- 
culated in the Fock basis, 



Z ss = E exp 



- > ( — n (n - 1) - pn 



(A9) 



At T = 0, the mean occupation number is then the inte- 
ger n that minimizes the exponent, 



U 2 



(A10) 



Following the continuous time path integral formalism 
yields the wrong result for the partition function, 



Z 'ss = E 



exp 



' ( + P n 



(All) 



and hence the wrong result of (n) — p/U. 

However, application of the discrete time path integral 
Eq. (A5| yields the correct value. Using a coherent state 



basis one finds 

L t = - log[{p t ip t | Pt+HPt+i)} 

{p t <fit\ H ss \p t+1 ip t+1 ) 



At- 



{pm I Pt+ivt+i) 

lipt + Pt^-s/m^e^^-^ 



(A12) 



UAt 



-ptPt+ie 



2i(<Pt+i—<p) 



- pAt^p^Pe^+^l 
where At = p>/N t . Using a saddleploint approximation, 
pt=p+6p t , (A13) 

we have 



Z = / VSpV<j) 



where 



exp - E £ o + C\ + £ 2 + £ 



£n = ( jp 2 -pp)Al. 



£1 = (Up-p) AtSpt, 



C 2 = - [1 - {2Up - p) At] x 



(Sp T+1 -8p T ) 2 



(A14) 



(A15) 



(A16) 



if 



+ p( 1 Pt+i - <Pt) 

i {Sp T+ i + Sp T ) (lpT+1 - tp T ) 

UAt (Sp 2 t + S P 2 t+1 ) . 



(A17) 
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We choose p = p/U, as in the continuous time inte- 
gral case, to minimize the zeroth-order action. Then in 
momentum space 



Z = J V5pV4> exp 

where the sum is over ui 
after the substitution, £q 



(A18) 



_ 2tt Nt — 1 2tt N t -1 it 

P 2 • ■ ' 2 ' ±±cic ' 

= and 



£2 = \ (6p u <p u ) GZ 1 (A19) 



where 



l G v ]l,l 



(1 - cos (uAt)) + (1 + cos (uiAt)) pAt 
2p/U 



[GJ 1 ] 1)2 = - [GJ 1 ] 21 = - sin (ojAt) (1 - //At) 
[GZ 1 ]^ = -2 (p/U) (1 - cos (wAf)) (1 - M At) . 

We invert G _1 to find the propagators 
{8p^5p v ) = 5 u - n n/U 



(A20) 



(A21) 



^LO. — rj 



u 
If, 



(1 + //At) - cos (wAt) (1 - //At) 
(1 - //At) (1 - cos (wAt)) 



Next we calculate 
OF 1 



dp 







9// 



9£| 
d// 



(A22) 



where we have neglected £ 



in 1 



Inserting the values 



Eqs. (A18)-(A19l, we find 



u~ 



E 

- £ 
" u 



-cos(cjAf) . 

(1-2/jAt) 
4p 



I At sin (wAt) cot (wAt/2) 



(l+/jAt)-cos(^At)(l- 
(1-jU.At) 



j£At) 



1 At 



1 



/3 ^ 2 1 - pAt 



p_ 

U 



1 



1 



2 1 - pAt 



As we take At — > this becomes 



(A23) 



§, in 



agreement with Eq. (A10). A subtle error still remains, 



namely that at zero temperature the occupation number 
must be an integer. To restore this constraint one would 
need to explicit sum over the topologically distinct "in- 
stanton" paths where ip has multiple windings. 
Appendix B: Explicit Calculations in the Discrete 
Time Path Integral 



We provide here a further explicit example of a calcu- 
lation in the discrete time step path integral formalism. 
For a more elementary example see Appendix [XJ The 
most important result in this appendix is the develop- 
ment of a formalism in which a discrete time calculation 
is expressed as the sum of a continuous time one and 
some easily calculated corrections. 



Our starting point is Eq. ( 21 ). Explicitly, the quadratic 
term is 



£ 



= \ ( w J***) g ^ ( 



\fp4>-k,- 



(Bl) 



of the propagators in Eq. (A22) into the derivative of where Q 1 is an inverse Green's function matrix, 



G. 1 



( 2(1- cos (a; At)) + 2 (1 + cos (wAt)) pV At 
+4 JAt cos (wAt) J2 d (1 - cos (k d a )) cos (A$ d ) 
+4i JAt sin (wAt) J2d sin (^rf a o) sin (A$d) 



2 sin (to At) (1 - pUAt) 
-4 JAt sin (uAt) J2 d i 1 ~ cos ( k dao)) cos (A$ d ) 
+4i J At cos (a; At) J2 d sin (k d a ) sin (A$ d ) 



-2 sin (w At) ( 1 - pUAt) \ 
+4 JAt sin (uiAt) J2d i 1 ~ cos (kd,a )) cos (A$ d ) 
—AiJAt cos (u)At) J^d sm (^d a o) sin (A^) 

2 (1 - cos (wAt)) (1 - pUAt) 
+4 JAt cos (wAt) (1 - cos {k d a )) cos (A$ d ) 
+4i JAt sin (a; At) sin (fc d a ) sin ( A$ d ) / 
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The propagators are obtained by inverting Q , and are given by 



(8pk,u<t>q,ri) 
{4>k,uj(t>q,ri) 



gtg) (q + k)6 u> . 
(a /27r) D 



_2 (1 - cos (wAt)) + [(fife + f 2fe ) cos {to At) + (£ lk -£ 2k )]At 

P 1 q t 1" t> (A<s) 

2(1- cos (uAt)) (1 - \ (£i k + £2k) At) + f fe 2 At 2 



(a /27r) u 
(ao/27r) D 



■ sin (wAt) 



1 - \ (£ lk + £ 2k ) At 



2(1- cos (wAt)) (1 - | (fife + f 2 fe) At) + f fc 2 At 2 

1 2 (1 - cos (a; At)) + [(f lfc + f 2fc ) cos (a; At) + (f 2fc - £ lk )} At 
4,p 2(1- cos (wAi)) (1 - I (fife + f 2 fe) At) + f fc 2 At 2 



0(A*) 
F0(A$) 



(B3) 



As explained below, the only important features of 
these functions are their ujAt —¥ structures and their 
values at |wAt| = ir. We illustrate this result by calcu- 
lating the average occupation number via 



(n) 



dT _ 1 1 dZ _ 

9/7 ~ /?y z"977 ~ 

£(£)7 dDk ? 



(B4) 



is performed over a circle of finite radius % ^ 1 < |w| < 
^ jVt 2 +1 ■ In terms of the integral over this contour 7, the 
summation over frequencies can be expressed as 

F(u) 



-y 



F{u>) 



1 
27 



dui- 



— i Res 



- 1 

F(w) 



1 



(B9) 



«o 



V /9£q 
D \ dp 



dc 



dp, 



dp 



As we have assigned p — jj (/i + 2J^ d cos (A$ d )), the 
derivatives are given by 



d_ _ d_ 
dp dp 



dp d 
dp dp 



(B5) 



We calculate this quantity at A3? = 0. 

The saddle point contribution comes from the constant 



dp 



= -pAt. 



The contribution from this term to Eq. (|B4| is 

D 



1 

w 



(s)7 iDl £ 



V /dCo 
dp 



D 



a 



(B6) 



(B7) 



The sum on the left hand side is over the frequencies 

u n = - ■■■ 2 f ' tne sum 011 tne ri S nt is over 

the poles ujp of F (w) inside the contour 7, and 7 is the 
complex circle defined by \w\ = ^f^r = St- The no- 
tation Res (X (w) , wp) refers to the residue of X (u>) at 

w = LUp. 

In a continuous time calculation, one takes the con- 
tour 7 to infinity. Assuming F (oj) is well behaved, the 
integral on the right-hand side of Eq. (B9| then vanishes. 



In our case we must explicitly include this term. To cal 



culate the contour integral, we take u = ^7 
X = . . . 2-7T. As At — > 0, the Bose factor is 







1 



-1 < x < k 

7T < X < 27T 



and so the integral of Eq. ( B9 1 becomes 



did- 



1 l At 



d X e ix F 



e lx , with 



(BIO) 



(Bll) 



The nontrivial part of the calculation comes from the 
erm involving £ 2 ""'> 

dC k ' Un 1 

— | — = - (1 - cos (wAt) + fife At cos (uiAt)) x 
9/t /i 



Spk,u;5p- 



+ sin (wAt) At6pk,w<f>-k,-u 
- (1 - cos (wAt)) 2Atp0 fc , t 



-fe,— oj • 



In the limit At — > 0, the poles of the functions in 



Eq. (B3 1 converge to finite values of U). Hence w^At <C 1, 



and the sum over ujp accesses only information about the 
low-energy structure of Eq. (B3) and (|B8|). 



Thus, the summation of functions of the propagators 
in Eq. ( B3 1 requires only the form of their low- frequency 



(B8) poles, at wAt <C 1 (marked with superscript p) and their 
values on the contour 7, at ui = j^.e tx (marked by super 
script o) . For our particular case, the summand Eq. (|B8 
is composed of 



We perform the summation over the frequencies w„ 
by taking a contour integral. The same trick is used in 
the continuous time approach, but the contour here is 
slightly different. As illustrated in Fig.[9j the integration 



dp 



fifcAt 
P 



p 

fc,CJ 



4p 



(B12) 



O (At) 2 x ((SpSp) , {8pct>) , (#)) , 
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Im(6>) 




Re(ai) 



FIG. 9: The contour 7 used to per form the summation over 



2tt JVt-1 2tt JV f -l 



„ — , , asinEq. (B9l. The contour is given 

by ui = ^e lx , x — . . . 27r. As one goes to the continuous 
time case with At — > 0, the radius of the contour goes to 
infinity. 



and 



fringe*) At(J^)° 



- (l - cos (7re lx )) 2Atp\^ L 
1 



(1 - cos (7re lx ) + £i fc Ai cos (ttc 1 *)) 



(B13) 



4p 



Thus we do not need the full structure given in 



Eq. (B3) byte rather just the pole and contour values 



given in Eq. ( |22[ ) 

We 
JLr k 



We now explicitly calculate the contribution of 
9 £^,w» ^ Q 'j'^g low-frequency behavior is 



0£i 



V 



'Ik 



dp 



with poles at up = H£k, so that 

„fc,u;\ p 



i Res 



/ 

while the contour value is 

/ OCn \ 



duj 



\ dp 



coth (/3£ fc /2) . 



V 1 A 

a? 2 



27T X, \ <9/i 



-ir l = 



1 y 

"2< 



Hence 



1 



ao 



/3V V2tt 

_ 1 
_ 2 



d D k 

(2vr) D 



1 _ coth ( /3ffc / 2 ) 



(B14) 



(B15) 



(B16) 



(B17) 



Combining this result with the zeroth-order contribu- 



tion in Eq. (B7| we find 



1 /-a?d D k 



(2tt) 

where p = + 2JD) /U 



l-^coth(/3£ fe /2) 



(B18) 
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